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Abstract 



We continue our consideration of a class of models describing the reversible dynamics 

2 ' of A^ Boolean variables, each with K inputs. We investigate in detail the behavior of the 
Hamming distance as well as of the distribution of orbit lengths as A^ and K are varied. We 

i-^ I present numerical evidence for a phase transition in the behavior of the Hamming distance 

Ch ' at a critical value Kc ~ 1.65 and also an analytic theory that yields the exact bounds on 

3 ■ l-5<Kc< 2. 

We also discuss the large oscillations that we observe in the Hamming distance for K < Kc 
CN . as a function of time as well as in the distribution of cycle lengths as a function of cycle 

length for moderate K both greater than and less than Kc- We propose that local structures, 
or subsets of spins whose dynamics are not fully coupled to the other spins in the system, 
^p . play a crucial role in generating these oscillations. The simplest of these structures are 

linear chains, called linkages, and rings, called circuits. We discuss the properties of the 
Q ' linkages in some detail, and sketch the properties of circuits. We argue that the observed 

^D ■ oscillation phenomena can be largely understood in terms of these local structures. 



H . Key words: Gene Regulatory networks; Random boolean networks; Time-reversible 

Boolean networks; Cellular automata; 



1 Introduction 

Kauffman nets [1] have been used to model the complex behavior of dynamical sys- 
tems ranging from gene regulatory systems [2-4], spin glasses [5,6], evolution [7], so- 
cial sciences [8-11], to the stock market [12]. A Kauffman net consists of A^ boolean 
variables or spins, each of which is either "-|-1" or "-1" at each time t = 0, 1, 2, . . . . 
We define aj to be the value of the i^^ spin at time t. The system evolves through a 
sequence of states S", S^, . . . , where 

i:' = ialai...,a%). (1) 

There are 2^ different possible S^'s. 

The Kauffman net evolves according to 

^l = F.{^'-') , (2) 

where Fi is a boolean function with K boolean arguments picked from among the N 
possible arguments in S*^^ and one boolean return. The Fi for i = 1,2, .. . , A^ are 
fixed during the evolution. The combination of them is called a realization. There are 

^2(2 ))^(^j different realizations for a net of A^ elements with K inputs [2]. 

After specifying a realization, one picks a starting state for the system; the network 
then evolves following equation (2). 

Reference [13] (which we shall call paper I) introduces a time-reversible [14-16] 
boolean network model, [17] with dynamics governed by the equation 

cr*+i=F,(SVr • (3) 

Each time-reversible network realization has a corresponding dissipative realization 
with same functions and connections. In the reversible model substates S*^^ and S* 
are both required to calculate S*"*"^ [18], so the state of the system is (S*^"*^, S*), and 
the state space has 2^^ points. 

Paper I studies how the distribution of cycle lengths generated by the time-reversible 
model scales with A^ for different values of K. This paper investigates the cycle length 
distribution in more detail and also investigates the behavior of the Hamming dis- 
tance [5,19,20] in this model. In the next section, we present both numerical evidence 
and analytic arguments that the behavior of the Hamming distance changes quali- 
tatively and a phase transition occurs in the model at a critical value K^ — 1.65. In 
section 3 we investigate how the realization average of the average number of cycles 
depends on the cycle length. For small K, cycle lengths divisible by integer powers of 
two and by three are much more likely than lengths that are products of larger prime 
numbers. Oscillations are also observed in the behavior of the Hamming distance as 
a function of time when K is small. We explain these observations in terms of "local 



structures," groups of spins whose dynamics do not depend on the state of any other 
spin in the system. In section 4, we discuss the two most important structures, link- 
ages and circuits. We calculate Hamming distances and the cycle-length distribution 
for linkages. Our calculations for these structures are then connected back to the 
numerical observations, both for K = 1 and for higher K. 



2 Hamming distance in reversible Boolean nets. 



In this paper we define the Hamming distance D{t) to be the number of spins in two 
substates that are different at time t. 

Consider two different time developments of a given realization of a system that 
are identical at t = and that have a single spin different at t = 1. We obtain 
two series of configurations, called two paths of the system, and calculate D{t), the 
Hamming distance between the two paths as a function of time t. In this subsection 
we consider the behavior of {D{t)), the Hamming distance averaged over different 
initial configurations and realizations, as a function of time t, for different values of 
the parameters K and N in our time-reversible model. We will find that when K 
is small, at long times the maximum value of the Hamming distance saturates at a 
value that is independent of A^ as A^ -^ oo, while when K is sufficiently large, the 
saturation value is proportional to iV as A^ ^ oo. These two behaviors are separated 
by a phase transition at Kc, a critical value of K. 

It is useful to define the normalized Hamming distance {d{t)) = {D{t)) /N, the av- 
erage over realizations of the fraction of variables having different values in the two 
paths at time t. Below the phase transition, one finds that as iV — > oo {d{t)) tends 
to zero for all t. while above K^ the maximum of {d{t)) has a non-zero limit. 



2. 1 Hamming distance for K = Q and K = N 



When i^ = 0, if one starts with two copies of the system that are identical at time 
t = and have exactly one spin different at time t = 1, the Hamming distance 
oscillates between and 1 at all future times. When K = N, flipping a single spin 
changes an input of every spin, so that all spins have probability 1/2 of being the same 
in the two paths at time t = 2. Thus, when K = N, for all times t > 2, (-D(t)) = A^/2, 
and {d{t)) = 1/2. Thus we see that the Hamming distance at long times (-Dqo) = 
limi_>oo -D(t), is independent of system size when K = and is proportional to the 
system size A^ when K is very large. 



2.2 Hamming distance for intermediate K 



To investigate intermediate K, we simulated reversible networks of 100, 1,000, 10,000 
and 100,000 boolean variables for different K values. For each case, 100 realizations 
were generated, and for each realization, one pair of paths with -D(O) = and -D(l) = 
1 was examined. Numerical results for A^ = 100,000 are presented in flgure 1. For 
all i^'s, the Hamming distance flrst increases, then reaches a plateau, and flnally 
oscillates about that plateau. For small K, large oscillations in {D{t)) are observed. 
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Fig. 1. Plot of {D{t)), the average over realizations of the Hartiming distance as a function 
of time t, for systems with N = 100, 000 and different values of K. The starting initial 
conditions for each realization were chosen so that -D(O) = 0, -D(l) = 1. For visual clarity, 
this figure is split into two parts with different vertical scales. The left subfigure includes 
only i^-values below or close to the critical value Kc — 1.65, while the right subfigure 
includes i^- values significantly above Kc- 



In particular, when K is small, {D{t)) is especially small when the time t is a multiple 
of 2™ for integer m (e.g., at t = 64, 128, and 256). For larger values of K these 
oscillations become much less pronounced. 

This figure shows that (D)^^^, the maximum value of the Hamming distance, exhibits 
a sharp jump as a function of K at Kc ^ 1.65. This jump arises because, for K > Kc, 
(-^)max i^ proportional to A^ for large A^, while for K < Kc, (-D)^^^ is independent 
of A^ as A^ — » oo. This behavior is shown in figure 2, which plots the maximum level 
of the normalized Hamming distance (d)^^ = (D)^^ IN as a function of K for 
different A^. Thus, dmax appears to act as an order parameter for a phase transition 
at Kc ^ 1.65. 



The phase transition that we observe at Kc — 1.65 is of a percolative nature [21-24]. 
When K > Kc, if one changes a few spins in the system, then information about the 
change spreads to more and more spins and eventually covers a nonzero fraction of all 
the spins in the system. When K < Kc, information remains localized within a few 
spins which tend to oscillate with periods which are powers of small prime integers. 
This phase transition is then connected with the percolation of information within 
the system. 

The percolative nature of the phase transition in the reversible model has similarities 
to that observed in the dissipative Kauffman model [5,25,26], but there are differences. 
Specifically, consider two paths that are the time histories of a system realization 
that is started with two initial conditions that are almost the same, but having a 
small number of spin- values differing or 'deviating' from one another. The original 




Fig. 2. Average of the peak of the normalized Hamming distance (d)jnax = i^)ma.x/^ ^^ 
a function of K for four decades of N values. For each set of N and K, we generated 100 
pairs of paths with initial separations D{0) = and -D(l) = 1. Each of them is derived 
from a different realization and initial condition. For each pair of paths we calculate .Dmax, 
the peak value of the Hamming distance, and then average this over realizations. The inset 
shows the behavior near criticality. We use these data to estimate that the critical value 
Kr is about 1.65. 



Kauffman model is dissipative, so a spin whose values are deviating at early stages 
of the history may stop deviating as time goes on. However, a reversible system must 
undergo cyclic behavior. If a spin deviates at one time, it will do so at an infinite 
number of later times. In a reversible system, deviation is a permanent condition. 

So we must ask how often a deviating spin triggers the deviation of others. Let G 
be the average number of additional spins whose deviations are directly caused by 
a lone deviation in the system. The value of G determines whether a few deviations 
will cause a chain reaction which will spread to the entire system. If G > 1, the chain 
reaction will occur and the Hamming distance will eventually saturate at a value 
proportional to A^. In the opposite case, G < 1, a single deviation triggers less than 
one additional deviation on average, and the effect saturates long before the entire 
system is affected. 



2.3 Bounds on the critical value Kr. 



Though we have not calculated exactly the branching ratio G for this chain reaction, 
we can prove that 1.5 < Kc < 2. We first prove that 1 < K^ < 2, and then refine the 
argument to show that 1.5 < Kc < 2. 
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Table 1 

The different classes of functions, F{ai,a2)^ for K = 2. The classes are defined by the 
symmetry operations of spin flip and interchange of the two spins. 



We define systems with fractional K to be composed of elements that have probability 
[K] + 1 — i^ of having [K] inputs and probability K — [K] of having [K] + 1 inputs, 
where [x\ is the largest integer less than or equal to x. (This then defines K as the 
average number of inputs, since {[K] + 1 — K) + {K — [K]) = 1 and [-ft']([-ft'] + 1 — 
K) + {[K] + 1){K - [K]) = K). The growth rate Gk is the average 



Gk = m + 1 - K)G[K] + {K- [K])G^K]+i 
Thus we need only calculate Gk for integer K. 



(4) 



2.3.1 



K,>1. 



We can see immediately that deviations will not percolate if K is less than or equal 
to one. When K = 0, every spin is independent of its inputs, so Gk=o=0. When 
K = 1, there are four different possibilities for Fi{a), namely 1, —1, a, and —a. In 
our model all functions occur with equal probability, so a deviation in the spin a 
triggers a deviation in the spin CTj to which it is an input with probability 1/2. Thus 
the growth factor for i^ = 1 is 



G 



K=l 



0.5 



(5) 



The growth factor for K between and 1 is just KGk=i, so for K less than or equal 
to one the system must be below its percolation threshold. 



2.3.2 



Kc<2. 



For any K, changing any one of the K arguments of a function -Fi(S) has a 50% 
chance of flipping its output value. Therefore, a deviation of a single spin will, at the 
very next step, cause an average of K/2 deviations. Later steps may cause further 



deviations, but considering only this initial step certainly bounds the growth factor: 

Gk > K/2 . (6) 

Thus, when K > 2 the growth factor Gk>2 > 1, so that the deviation will spread to 
a nonzero fraction of the entire system. 

At first sight, it appears that i^ = 2 is the marginal, or critical case. That is wrong. 
Information will also percolate through the entire system when K = 2.ln the previous 
paragraph, we considered the further deviations induced in the very next step after 
the deviation of a given spin. These effects were enough to ensure that at least 
K = 2 is marginal. However, a deviation can have a delayed further effect, producing 
additional deviations beyond the one we have already calculated. These additional 
deviations then permit information to percolate through the entire system for K = 2. 

To see this result, we need only establish one case with a further, delayed, deviation. 
Notice Table 1. Consider a situation in which spin number 3 is determined by the 
"or" function shown on the last line in the table. Assume that input spin 1 has just 
deviated and input spin 2 has the value +1 on both paths. Now spin 3 does not 
deviate at this step. However, if at a later stage spin 2 takes on the value —1, that 
will induce a delayed deviation of spin 3. Since situations like these occur with non- 
zero probability, and increase the value of Gk=2 beyond the value shown in equation 
6, the case i^ = 2 is above the percolation threshold. 

2.3.3 Kc > 1.5. 

Since we know 1 < Kc < 2, we need only consider values of K between 1 and 2. We 
have already calculated Gk=i and a lower bound to Gk=2', now we wish to consider 
Gk=2 in more detail. 

Let us estimate the average number of additional deviations produced by one deviated 
spin. The growth rate Gk=2 is 2N times the probability that a spin a is induced to 
deviate when one of its two inputs has deviated, which is in turn the sum of products 
of the likelihood of a function choice times the probability that the function gives dif- 
ferent outputs with one deviated input. Of the sixteen two-input functions -F(o"i, (T2), 
two functions (±1) depend on zero inputs, four functions (iai, ±0-2) depend on one 
input (they change output if one input is changed but not the other one), and two 
functions (±o"i(T2) depend on two inputs (they change if either input is changed). The 
remaining eight functions are like the "or" function; they depend upon one argument 
only for a particular value of the other argument. (See Table 1.) Therefore, 

2x0+4x1+2 x2+8xg 

Gk=2 — TT. , 

lo 

where q is the unknown effective number of inputs of the functions like "or". We 
overestimate q by saying that these functions all behave exactly as if they had two 



inputs. This estimate then gives us 

Gk=2 <\. (7) 

Putting together our results from equations 5, 4 and 7, we find that 

G'< (2-ir)0.5 + (ir-l)1.5 = -0.5 + ir. (8) 

This overestimate gives us a marginal growth rate at K = 1.5, which yields a lower 
bound to the critical value Kc > 1.5. Our best numerical evidence is computed from 
the data shown in figure 2 and gives Kc = 1.65 ± 0.10. This result is, of course, in 
agreement with the exact bounds. 



3 Oscillations in orbit period distributions 



In this section we present our numerical results for cycle length distributions, and 
give a qualitative explanation of what we are seeing. When K is small but larger 
than Kci the number of cycles does not depend smoothly on length. Instead, cycles 
whose lengths are divisible by small integer powers of two and by three predominate. 
As K is increased, the cycle length distribution becomes smoother. We relate this 
trend to the evolution in the connectivity of the network as K is increased. When 
K is small, we can have small clusters of spins that are only weakly coupled to the 
other spins in the system. These local structures are shown to have a characteristic 
effect of yielding cycle lengths divisible by small factors such as two and three. As K 
increases, all the spins tend to be in a single large cluster with complex interactions, 
and local structures become much less likely. This leads to a much smoother cycle 
length distribution. 



3. 1 Simulation results for cycle length distributions 



Figure 3 shows the distribution of cycle lengths for intermediate K values and A^ = 10. 
The distribution oscillates strongly when K is small. As K is increased, the oscil- 
lations diminish and the distribution eventually converges to the random hopping 
behavior discussed in Paper I [13]. A complex oscillation pattern in the distribution 
is evident, particularly for the cases K = 2 and K = 3. Several different branches 
are observed, which reflect periods which are or are not divisible by small powers of 
two and of three. The oscillations at moderate K in the reversible model are signif- 
icantly more pronounced than those observed in dissipative Kauffman models [27]. 
The amplitude of the oscillation diminishes as K increases, and as K approaches A^, 
the distribution approaches the two-branch structure (due to an even-odd oscillation) 
calculated in the K = N case in Paper I. 

To see the oscillatory structure in more detail, we show in Fig. 4 an enlarged plot 
which shows the number of cycles for the region extending from / = 360 to / = 408. 
The most likely cycles fall into three main branches, which in decreasing order of 
probability satisfy / = mod 4, / = 2 mod 4, and / = 1 mod 2. Each main branch is 
further split, with cycle- lengths that are multiples of three being more likely. Thus, 
the cycle lengths that are divisible by 24 occur with the highest likelihood. 

This oscillation is not a small effect. Figure 4 demonstrates that cycle lengths divisible 
by 24 are more than two orders of magnitude more common than any odd cycle length 
in this range of /. Thus we see a marked tendency for orbit periods to contain factors 
of 2" for small integer n and of 3. 



10 



10 



I I ^1 llllll — ^ I |l llllll TTTTTTg 



10 



I I I llffll • I 
■ be 



I I I llllll I TTTTTB 










I 10 

fe 10 
> 

< 



(C)K=4 ♦ 

Ii Ii I I itiiii 



10 



10 



3 4^" 12 3 4 

10 10 10 10 10 10 



Cycle Length I 



Fig. 3. Log-log plots of numbers of cycles per realization, averaged over realizations, as a 
function of cycle length / for A^ = 10 and for K = 2, 3, 4, and 5. For each K, all cycles 
are enumerated from each of 25600 realizations. Note the different branches arising from 
oscillations in these functions. The different point styles distinguish the different divisors 
of the /-values. They are: I = mod 4 (circles), I = 2 mod 4 (squares) and I = 1 mod 2 
(diamonds). Filled symbols have I = mod 3, while open symbols denote I not divisible by 
three. 



3.2 Characterizing the oscillations in the cycle length distribution 



We now look for a quantitative measure of the amount of osciUation in cycle length, 
so that we can see how it changes as a function of K. More specifically, we wish to 
ask whether we can understand how the oscillations go away for large K. 

For small K, the even cycles predominate. It is very likely indeed, that among the 
many independent localized stuctures in oscillation at least one will give a period 
two oscillation, which will then be observed in the overall period of the system. To 
show this, we plot in Figure 5 the ratio of the number of even cycles to the number of 
odd ones for two different system sizes, A^ = 10 and N = 7. For large K, the global 
behavior forces the number of even cycles to be greater than the number of odd ones. 



11 



10' 



10 



10" 



u 



10 



1 10- 
lO' 
10" 
lO" 



• 0, 12 mod 24 

o4, 8, 16, 20 mod 24 

-■6, 18 mod 24 

-a 2, 10, 14, 22 mod 24 

-♦3,9, 15, 21 mod 24 

-«1,5, 7. 11, 13. 17, 19. 23 mod 24 



380 384 388 

Cycle Length Z 



Fig. 4. Enlarged plot on a semilog scale of the realization average of the number of cycles 
of length / versus / for A^ = 10 and for K = 2. The oscillation structure is highlighted by 
showing only the region extending from / = 360 to / = 408. The main structure of the figure 
consists of three branches, which in decreasing order of probability satisfy: / = mod 4 
(circles), / = 2 mod 4 (squares) and / = 1 mod 2 (diamonds). Filled symbols are used for 
lengths evenly divisible by three and open symbols for those that are not. The plotting 
symbols are the same as in Fig. 3. 
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Fig. 5. The ratio of the number of even cycles to that of odd ones plotted against K for 
N = 7 and N = 10. For each K and system size, all cycles are enumerated in each of 25, 600 
realizations. For small K values, the ratio increases rapidly as the system gets larger, while 
for large K, the ratio stays at a level above unity, which depends much more weakly on 
system size. 
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Fig. 6. Q{r) as a function of r. Top: even cycles. Bottom: odd cycles. High Q values mean 
high likelihood of observing a cycle length with certain residues. Numerical data were 
obtained by complete enumeration of 25,600 realizations for each K. 



but the effect depends very weakly on A^. For smaller values of K the ratio becomes 
much larger, and has substantial A^-dependence. 

To parametrize the oscillation structures for different K values, we first separate the 
domain of possible cycle lengths into segments of length P. Based on the argument 
about the appearance of factors of small powers of 2 and 3 in the system's periods, 
and also the observation in Fig 4, we choose P to be 24. We also require the ith 
segment to start at 24i and thus the 24 elements of the segment are 24i + r, r = 
0,1,... , 23. Paper I discusses how correlated fiuctuations in the number of mirrors 
could produce an even-odd oscillation in the number of cycles. In this paper, the even- 
odd oscillations (exhibited in figure 4 of Paper I) are just a distraction. To eliminate 
the even-odd effect due to the correlated fiuctuations in the number of mirrors, we 
define qi{r) to be the proportion of cycles of period 2Ai + r among all the even cycles 
in the segment when r is even, and among all the odd ones when r is odd. We then 
average qi{r) over segments of which the number of both even cycles and odd cycles 
are observed more than 100 times in our simulation (complete enumeration of 25,600 
realizations), and thereby obtain the functions Q(r). 

The functions Q{r) are plotted separately for even and odd values of r in figure 6. 
For K = 2 we can see very strong oscillations as a function of r in both plots. As 
expected, factors of three and of four are favored (this plotting method does not show 
the oscillations with period two). As K is increased, the r-dependence becomes much 
weaker. 

To characterize the strength of the oscillation for a given K using a single parameter. 
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Fig. 7. Our Gini coefficients measure the size of oscillations in the number of cycles as a 
function of length. They are obtained by complete enumeration of 25,600 realizations for 
each K and plotted here, with even and odd cycles separated. When K = 1, the only odd 
cycles observed are of length 3 and thus G = 1 for the odd cycles. 



we order r = 0, 1, 2, ■ ■ ■ , P — 1 in the increasing order with respect to Q{r) and index 
them as ri,r2, ...,rp. The Gini coefficient [28], a measure of the inequahty between 
the P values of Q{r), is defined as 



^ = 1 



(p-i)Ef=ig(r,; 



(9) 



When all the Q{r) are equal, G takes on its minimum value, G = 0; when only Q{rp) 
is nonzero, so that the distribution of Q{r) is extremely inequitable, then G = 1. 

The Gini coefficients as a function of K with A^ = 10 and P = 24 are shown in 
figure 7. Notice that there are separate curves for even and odd values of r. We see 
that the Gini coefficient is large for small K and decreases rapidly as K gets larger, 
showing values indistinguishable from zero for K above 5. 



3.3 Local structures 



In this section, we discuss why cycle lengths are so likely to be evenly divisible by 
small powers of two and by three. When K > K^, the influence of a single spin 
flip percolates and reaches a flnite fraction of the entire spins in the lattice. We can 
then imagine that a nonzero fraction of all the spins in the system fall into a main 
connected cluster, one with many complex linkages among the different spins in the 
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cluster. We depict this in Figure 8, in which the 'C denotes the connected cluster. It 
is further indicated by the surrounding curve. 



3.3.1 Isolated structures 

But K is not infinite. There will be spins in the system which are only weakly coupled 
in to the rest. For example, there will be spins which are simply uncoupled to any 
others. One such spin is indicated by the 'Jq' in the figure. This spin has a function 
Fj(S) which is simply ±1, independent of all the spins. Thus, this coupling has an 
effective i^-value of zero. It is possible that this isolated spin does not affect any 
other spins in the system. It can also happen that the isolated spin can affect one 
or several spins in the cluster. We show this with the situation labeled as 'Ji'. The 
effect of the two situations is rather similar. 

There are other examples of weakly coupled spins. But the effect of the isolated spin 
can be calculated in full detail, and will serve to illustrate the qualitative nature of 
localized-structure-effects. 

First of all notice that for large N there will be many such isolated spins. For a given 
value of K the probability that a spin is assigned the function —1 or +1 is 2*^^~^ \ 
so that the average number of such spins will be 



P 



)(l-2-)^ . 



This number must be large for sufficiently large A^. A spin with input function —1 is 
in a period-4 cycle, while a spin with input function +1 is in a cycle of period 1 for 
half its possible initial conditions and in a cycle of period 2 for the other half. Thus, 
for large N there will on average be many isolated spins in cycles of length two or 
four. It follows that it will be overwhelmingly likely that the cycle length of the entire 
system will be divisible by two. 

Our A^'s are not very large. But A^ = 10 is large enough so that isolated spins are 
likely, and they play a substantial role in making our cycle-lengths mostly divisible 
by several factors of two. As figure 5 demonstrates, when K is small, the ratio of even 
to odd cycles increases rapidly as A^ is increased, consistent with the idea that local 
structures have a substantial effect at small K. 

Factors of three are only slightly harder to obtain. Imagine, for example, a spin 
coupled only to itself via 

F,(S) = ±ai . 

One initial condition for such a spin results in a cycle length of one, while the other 
three result in a cycle length of three. These and other structures can then give us 
our observed factors of three. 

Other isolated structures can be formed that have a substantial influence on system 
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C: Connected Cluster I: Isolated Spins 



L: Linked Spins 



D: Dangled Spins 



Fig. 8. This figure shows the simplest kinds of clusters. The arrows show the flow of in- 
formation. The cluster inside the circle labeled 'C is connected. That means it contains 
many spins intimately linked by complex interactions. The clusters labeled '/q', '/i' and 
'I2' are single spins, isolated from the influence of any others. As we have depicted them, 
the '/q' conveys no information to the connected cluster while '/i' brings information from 
the isolated cluster to the connected cluster. We also depict the linkages 'Lq' and 'Li', 
slightly more complex structures which also undergo their own independent motion. The 
D's indicate dangling structures. 



properties. For example, one can have a structure like the linkage 'Lq' shown in Fig. 
8. A linkage is a sequence of spins linked together by coupling of the form ±0". The 
arrow indicates that the value of the spin on at the tail changes the values of the 
spin at the head. A linkage like 'Lq' has no effect on the cluster. One can also have 
a linkage like 'Li' which directly influences the cluster. We shall discuss the effect of 
linkages in more detail in the next chapter. For now we merely point out that these, 
and other, more complex structures have a simple effect upon the overall periods of 
the system. Each small isolated structure has its own short period. The motion for 
the entire system must be a multiple of the periods for each of the isolated structures. 
Since the small structures all tend to have periods two and three and four, it is no 
wonder that these factors appear very often in the observed periods of Fig. 3. 



3.3.2 Dangling structures 

It is also possible to have isolated structures which are influenced by the connected 
cluster but do not influence it themselves. For example, consider a spin, at, whose 
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input function -Fj(S) depends only on spins that are not causally related to (Tj. Then 
the extra spin 'dangling on' the cluster produces a simple effect upon the periods of 
the entire system. To deal with this situation, we prove a theorem. 
Theorem 1. Consider two sequences F* and a*. We assume that both of them are 
sequences of +l's and —I's, that F* is periodic with period Tp, and that the a's 
depend upon the F 's according to 



(jt+^ = F' ■ a'-\ (10) 

Then a is periodic with period T„, where T„ can only be Tp, 2Tp, or 4Tp. 

We prove this by showing that (A) T^ divides 4Tp evenly, and (B) Tp divides T^ 
evenly. To prove (A), note that 

a*+2T^=n(t) a' , (11) 

where 

■Q/^\ ^ pt+2Tp~lpt+2Tp--i pt+-ipt+l Q2) 

When Tp is odd, n(t) is a product over a full period of the F's. Therefore, it is 
independent of time. If it is +1, the sequence of a's repeats after 2Tp. If it is -1, then 
a^ repeats after 4Tp. IfTp is even, then 



n(t) = [C^ii'F*+'^^-2-+i 



2 

= 1 . 



Therefore, when Tp is even, 0"*+^^^ = or*. Thus, for both even and odd Tp, a is 
periodic with a period T„ that divides ATp evenly. 

Now we prove (B), that Tp divides Trj evenly. For any t, we have 

a^-+* = n'(t) CT* with 

n'(t) = pT„+t-i . . . pt+zpt+i _ 

The periodicity of a implies that 11' (t) = 1 for all t. Note that, for all t, 

n'(t + 2) = F^''+*+i F*+i n'(t) . 

Thus, _F'^<^+* = _F* for all t, and Tp must divide T^r evenly. 
Q.E.D. 
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Theorem 1 enables us to see additional mechanisms for the multiplication of cycle 
lengths by factors of two and four. Look at the structures dangling from the end of 
a connected cluster as also shown in Fig. 8. A dangling cluster consisting of a single 
spin can multiply the periods produced by the connected cluster by a factor of one, 
two, or four. The two-spin dangle can, at most, produce a lengthening by a factor 
of eight. Thus, these dangles produce additional mechanisms for the appearance of 
factors of two in the period of the entire system. 
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4 Some local structures 



In the preceding two sections we saw that small K reversible Boolean nets exhibit 
large oscillations, both in the behavior of the Hamming distance as a function of time 
and the orbit period distribution as a function of period /. We explained the latter 
as an effect of small localized structures. In this section, we discuss the localized 
structures in more detail. Our purpose is both to explain the oscillations exhibited 
above, as well as to explore some of the structures' statistical properties. 



4.1 The structures and their coupling 



Two kinds of elements, depicted in figure 9, serve as a basis for a description of the 
behavior dX K = 1. These same elements, called linkages and circuits, may also be 
expected to dominate the behavior at small K. Others have discussed the effect of 
structures like these for the dissipative Kauffman model[26, 29-31]. The behavior is, 
of course, different in the reversible case. 



4-1.1 The structures 

A linkage is L spins, Uj, j = 1, . . . , L, coupled together in a linear array. The first 
spin has an input function, Fi, that is either +1 or —1. Information is transferred 
from one spin to the next, for all j bigger than one, with a coupling function 

F,(S) = ±a,_i for j = 2,3,...,L. (13) 

A circuit of length L is like a linkage except that the last spin is coupled to the first 
with 

Fi(S) = ±aL. (14) 

These basic structures may also be tied together, as for example, in the tree and the 
tadpole shown in figure 10. In this section, we shall analyze the basic structures. We 
shall, however, not explore the effects of the multiplicity of structures which can be 
present at one time, nor of the different ways these structures can be coupled. 



4.1.2 Local structures at K = 1. 

For K = 1, since two of the four Boolean functions of one variable do not depend 
on the input, half of all spins have dynamics that do not depend on the state of any 
other spin. Therefore, typical local structures in networks with K = 1 tend to contain 
only a few spins. 
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Linkage of Length 4 




Circuit of 
Length 6 



Fig. 9. Two kinds of basic structures. Each circle represents a spin; the spin at the end of 
an arrow depends on the spin from which the arrow starts. 





r~\ Tadpole 



Tree 



o 



Fig. 10. Possible ways in which the basic structures may be tied together. Both the tree 
and the tadpole are possible structures at K=l. However, if the arrows were reversed, these 
structures would only become possible at higher K. 



For a given K, the numbers of linkages and circuits depend only on the placement 
of the couplings and on the function choices, and do not depend on whether the dy- 
namical equation is dissipative or reversible. Flyvbjerg and Kjaer [30] have calculated 
the number of circuits of length L exactly in the context of the dissipative Kauffman 
model. However, they did not calculate the number of linkages, doubtless because 
linkages of all lengths have trivial dynamics in the dissipative Kauffman model. 

The number of spins in linkages is of order A^, while the number of spins in circuits is 
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of order unity. This scaling follows readily via the following argument. When K = 1 
each spin has one input chosen from N possibilities, and the spin depends on the 
input only if its function is icr, which occurs with probability 1/2. The probability 
that a spin depends on any other spin is of order one, while the probability that 
it depends upon itself is of order 1/N. Consequently, there are of order A^ spins 
linked to other spins, but only of order one spin linked to itself. An extension of this 
argument implies that there are of order A^ linkages, but only of order one circuit 
in the entire system. Tree-like structures can be formed from linkages alone, and 
tadpoles from combinations of linkages and circuits. The former have a behavior 
which is qualitatively similar to a single linkage, the latter have properties which are 
closer to that of circuits. 



4.2 Theorems for the analysis of linkages and circuits. 



There are three theorems which greatly facilitate our work with linkages and circuits. 

We have already stated one of these as our theorem 1 in Sec. 3. This theorem tells 
us how the period of a base structure is affected by a spin dangling on the end. The 
period of the base is multiplied by one or two if the base period is even and one, two, 
or four if it is odd. By direct calculation, one sees that a linkage of length one will 
have a period of one or two or four. Hence, all linkages have periods that are integer 
powers of two. 

The other two theorems apply more specifically to linkages and circuits. They are: 
Theorem 2 (Equivalent Structure Theorem). For a linkage or a circuit, switch- 
ing the function assignment Fj{a) from a to —a does not affect the dynamics of the 
structure in the sense that we can always find a one-to-one map from the state-space 
of the local structure to itself, relating a cycle before switching the function assignment 
to a cycle with the same cycle length after the switch. 

Proof: For each Fj{a) which has the value —a, change the variable crj_i into —aj^i. 
Q.E.D. 

Note: This change leaves all Hamming distances and cycle lengths unchanged. Hence, 
for the examination of these quantities we need only consider the case in which the 
linkage functions are -Fj(S) = (Tj-i and not the cases in which -Fj(S) = —aj^i. 

Theorem 3 is a superposition principle that applies to systems with K = 1. 
Theorem 3 (Superposition Principle). Consider a network with K = 1. The 
equations of motion are of the form 

af' = F,(S*)af 1 , (15) 

Let 
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fj = 1 if FjiTj^) = +1 or a\ for some i{j) 

fj = —l if Fj(S*) = — 1 or — o-\ for some i{j) . (16) 

Let (T* and /i* both he solutions to the equations of motion for the same realization of 
the K = 1 couplings. Then 



P] = fj ^/i5 (17) 



is also a solution. 



Proof: This multiplication rule can he immediately verified on each hond in equa- 
tion 15. Therefore it must he true glohally. 

Note that the change of variables from a to (In/ a) /{in) converts this multiphcative 
superposition into the kind of additive superposition principle which arises for all 
linear equations. 



4-3 Hamming distance for linkages and circuits. 



One can make considerable progress in the Hamming distance problem for long link- 
ages and circuits. The Equivalent Structure Theorem implies that for each L we need 
to study only one circuit problem and two linkage problems. The latter differ by their 
first coupling, which can take on the two values Fi = ±1. 

First we show that the Hamming distance can be obtained by solving a single linkage 
or a circuit problem. Consider a linkage or circuit with a solution S* = {cr*}, and let 
Ai^ = {/U*} be another solution with initial conditions that are identical to those for 
S* at t = and have, at t = 1, a single spin different, the spin with index J. Then 
defining 

P] = ^>5 , (18) 

all the values of p^ and pj are unity except for p^j = — 1 . The variable p* is minus one 
whenever the two paths differ, so the Hamming distance between Ai and S at time 
t is 

D{t)=j:la-p-)- (19) 

The subsequent history of p is calculated starting from the initial data and the equa- 
tions of motion 

P'j^' = pUpT' fo^ t = 2,3,...,L. (20) 
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Fig. 11. Cascading Effect of a Single Spin Flip. Here, /9*- is plotted for a situation in which 



the flip occurs in a very large linkage or circuit. The grey squares denote p 
p=-l. 



1; the black, 



For the linkage, all the p* are one if j < J and (—1)* if j = J . The picture just stops 
short aX i = L; nothing is meaningful to the right of this value of i. For the circuit, 
we apply the periodic boundary conditions p*^j^ = p*. 

The result for a very large linkage or circuit is shown in Figure 11. This figure shows 
the familiar Sierpinski gasket [32,33], rotated ninety degrees from the conventional 
representation. In fact, equation 20 is equivalent to rule 90 in Wolfram's cellular 
automaton scheme, which is known to generate a Sierpinski gasket [34]. 

Figure 12 shows a plot of Hamming distance D{t) versus iteration number t for the 
gasket. The distance is unity at times 2" and its average in the binary decade [2", 2"^^] 
is (|j . Thus, this average increases as t*^'°S2 3-i) ~ t"-^^ [35]. The maximum value 
of D{t) in the range [2"~^,2"] is the Fibonacci number Fn, where Fi = 1, F2 = 2, 



and Fj = Fj_i + Fj_2, which scales as t ""'V ^ J ^ t°-^^. Both the average of D{t) 
and the maximum value of D{t) in each binary decade can be obtained analytically 
by exploiting the fractality and symmetry properties of the Sierpinski gasket. 

To calculate the Hamming distance in a linkage of finite length L, started by fiipping 
a spin J sites into the linkage, note that no information travels to sites (Tj with i < J, 
and that the dynamics are cut off at the last site i = L. If the distance from J to the 
right hand boundary is between 2" and 2""*"^, the motion repeats with a period 2'^"'"^. 
We shall see below how this information may be converted into statements about the 
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Fig. 12. Numerically calculated Hamming distance versus time for a linkage with Sierpinski 
gasket dynamics (solid line) . The average of the Hamming distance in each binary decade 
grows as t(^°2 3) ^ ^0.59 ((Jotted line), while the maximum value of the Hamming distance 



in each binary decade grows as t °^- 



l + V^ 



(dashed line). 



distribution of periods for the linkage. 

The Hamming distance behavior for the circuits can be understood by seeing the 
circuit picture as one in which the gasket picture is folded onto itself an infinite 
number of times, with the folding distance being, L, the size of the circuit. Because 
two superposed black dots cancel and make a gray one, the folding can produce 
complex patterns. Nonetheless, we know that the Hamming distance is unity at times 
which are of the form 2", for integer n, since at these times there is but one black dot. 
Furthermore, the typical or average Hamming distance must grow no more rapidly 
than t"'^^, since this growth is bounded above by the number of black dots in the 
non-overlapping gasket. We shall come back to this fold-over picture in order to assess 
periods of circuit motion. 



4-4 Cycle length distributions for linkages. 



Table 2 shows the distribution of cycle-lengths for the first few linkages. As expected, 
the linkages all have periods which are powers of two. We now calculate the distri- 
bution of cycle-lengths. 

First we show that the orbit period of a linkage is determined entirely by the root spin, 
defined to be the spin aj with the smallest value of J that has an initial condition 
other than (+1,+1). It is useful to define the ( — h) primitive sequence ^^^^Sj^ for 
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number of cycles with period 


1=1 


/=2 


1=4 


/=8 


1=16 


1 




2 










2 




2 




3 






3 




2 




3 


6 




4 




2 




3 


30 




5 




2 




3 


30 


48 


6 




2 




3 


30 


240 


7 




2 




3 


30 


1008 


8 




2 




3 


30 


4080 



Table 2 

Numbers of orbits of different periods for linkages of lengths L from 1 to 8 with Fi = +1. 



L 


Cycle Length Vs (Multiplicities) 


1 


1(1) 3(1) 


2 


1(1) 3(1) 6(2) 


3 


1(1) 3(1) 5(3) 15(3) 


4 


1(1) 3(1) 6(2) 12(20) 


5 


1(1) 3(1) 17(15) 51(15) 


6 


1(1) 3(1) 5(3) 6(2) 10(24) 15(3) 30(126) 


7 


1(1) 3(1) 7(9) 9(28) 21(9) 63(252) 


8 


1(1) 3(1) 6(2) 12(20) 24(2720) 



Table 3 

Complete enumeration of all possible orbits for circuits of lengths L < 8. The numbers 
in the body of the table are the possible cycle lengths and, in parentheses, the number of 
different cycles of that length. 



a linkage of length L to be the values {o"*} for i = 1,....L, t = 2,3,.... that one 
gets starting a linkage of length L with Fi = 1 from the initial condition with all 
(T° = al = 1 for all i except for a^ = —1. The infinite primitive sequence •^"+^5* is the 
L -^ oo limit of the ^ ^^<Sj^- The (H — ) primitive sequence is obtained by starting a 
linkage with Fi = 1 with the initial conditions a^ and aj = 1 except for al = —1. For 
all t > 1, the (H — ) primitive sequence is obtained from the ( — h) primitive sequence 
by changing the origin of time by one. The primitive sequence ^ ^5* can be obtained 
by superposing ^^ ^iS* and *■ ^^5*. The primitive sequences ^~~^^Sj^, *■ ^^Sj^ and ^^~^^Sj^ 
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Fig. 13. Primitive sequences ^"^ ^5r^, ^ Sj^^ and ^^^^^Sj^ for a linkage of length L = 20. 



for a linkage of length L = 20 are shown in figure 13. We use S^ to denote any one 
of these three primitive sequences. 



As demonstrated in fig. 13, the l^^ spin (/ > 1) in the 
with the period 



primitive sequence cycles 



(21) 



where [x] denotes the greatest integer less than or equal to x. Obviously the period 
of the l*^ spin in the (H — ) primitive sequence is also p{l), and it is straightforward 

to verify that in the ( ) primitive sequence the l^^ spin cycles with the period p{l) 

if /> 1. 

Now we show that changing the initial conditions of spins other than the root spin 
does not change the period p{l). We do this by noting that the time development 
from any initial condition can be written as the sum of primitive sequences with 
different spatial origins. The form of p{l) guarantees that any superposed structure 
has a period that is either equal to p{l) or else divides p{l) evenly. Since the time 
development of the j*'* spin can be written as a superposition a* = A*i?*, where 
^t+p(0 _ j^t g^j^^ Qt+p{i) _ j^t^ clearly T^, the period of the a*, must either be equal 
to p{l) or else divide it evenly. 

However, we still need to show that T^- is not smaller than p{l) (that this is not 
trivial can be seen by noting that superposing the two period-four sequences (+ + 
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h H . . . ) and ( \- -\ 1- + ...) yields a period-one sequence). First 

note that the period cannot decrease if Tg, the period of the S*, is less than Ta, 
the period of the A*. This follows because if Tb < Ta, then 5*+'^'4/2 _ j^t^ ^^^^ 
^t+TA/2 ^ ^t+TA/2Qt+TA/2 ^ j^t+TA/2-Qt ^ j^t^t Therefore, we need only consider the 
case in which the periods Ta and Tb are equal. 

We use the Sierpinski gasket form of the primitive sequences to show that superposing 
sequences A and B does not decrease the period when / > 2. The spins in a primitive 
sequence which have equal periods have time evolutions that map out triangular 
shapes in the i — t plane. For example, in the (H — ) primitive sequence the values 
a\ in which i = 2" + j and t = 3 ■ 2" — j with < j < 2" for any integer n have 
a\ = —1. Superposing a sequence with a larger value of i in the range [2" + 1,2^"] 
cannot eliminate the —1 for the smallest value of i in this range. Since a\ = 1 for 
all t with 3 ■ 2" + 2 < t < 5-2" when i obeys 2" < i < 2^", combining any of 
these sequences cannot reduce the period. A similar argument applies to the ( — h) 
primitive sequence using the lower edge of triangles in the gasket. 

Thus we have shown that if a linkage has Fi = 1 and the root spin has an initial 
condition other than (+1,+1), the cycle length of the /*^ spin with / > 2 is p{l) = 

2[i°g2{'-i)+2]^ Since the spin a^ in the ( ) primitive sequence has the time evolution 

(+1,+1,— 1,— 1,...), which is the same time development of the root spin al in a 
linkage with Fi = —1, it follows immediately that the cycle length of the /*^ spin in 
a linkage with Fi = — 1 is 

P(l) = 2[i°S2(0]+2 _ (22) 



Now we are in a position to find the distribution of orbit lengths for different linkages. 
Let Nf-^=_i{L,1) and Nf-^^=i{L,1) be the number of orbits of length / in a linkage of 
length L with input functions —1 and +1, respectively. Calculating Nf^=_i{L,1) is 
simple. We know that all the orbits in this linkage have length P{L), and that there 
are 2^^ points in the phase space altogether, so the orbit length distribution is 

22L 

NF,=.,{L,l) = Si,p^L)j^y (23) 

To obtain the distribution of cycle lengths for a linkage of length L with Fi = 1, 
we must consider initial conditions in which o"° = (jJ = +1, where the behavior is 
identical to that of a shorter linkage. The three other initial conditions for o"i yield 
the orbit period p{L), so Nf^=i{L,1) satisfies the recursion relation 

Nf,=i{L,1) = 3^-—^6i,p^L) + Nf,=i{L - 1,1) . (24) 

p{L) 

The solution to this recursion relation has A^i?^=i(l,/) = 26i^i + 61^2, Nfj^=i{2,1) = 
25/1 + 6,2 + 3(5/4, and, for L > 2, 
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<P>=(5.2±0.1)ln(AO+(0.0±0.7) 



100 1000 

Size of the System N 



100000 



Fig. 14. Average cycle length (P) plotted against the size of the system N for the reversible 
model with K = 1. The average is over 100 randomly chosen realizations, with 100 cy- 
cles randomly generated for each realization. The solid line is a fit to the numerical data 
assuming a logarithmic dependence: (P) = (5.2 zb 0.1)lnA^ + (0.0 =b 0.7). 



/,(2[l°g2(i-l)l_l) 
nr (T r\ Z_ ^ /.L-2[l°S2(i-l)l -.W 



2[log2(L-l)] 

[log2{i- 



og2{-f^-l)]-l 42"-l 



(4^" - l)<5,,2"+2 . 



(25) 



n=0 



This result agrees with our siinulational data shown in Table 2. 



We can connect these results to the way the typical cycle length scales with system 
size N for K = 1. Since the realization average of the number of circuits per realization 
does not change with A^, we expect the linkages to dominate the scaling, with the 
circuits only giving a constant multiplicative factor in the typical cycle length. For 
K = 1, the probability of finding a linkage of length L decreases exponentially as L 
increases. Therefore, Lmax{N), the longest linkage in a system of A^ spins, obeys 



Lmax{N) OC l0g{N) . 

We just argued that the typical cycle length of a linkage of length L is 

P(L) oc2['°S2i] . 



(26) 



(27) 



Eqs. (26) and (27) together imply that the typical length of the cycles for i^ = 1 
should increase logarithmically with N. Our simulation results, shown in figure 14, 
are consistent with this scaling. 
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(A) Spin; (B) Spiny 
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Fig. 15. Dynamical development of deviations in a linkage (A) and in a circuit (B). 



4-5 Dynamics of circuits 



Table 3 gives the distribution of cycle lengths for circuits for the lowest values of L. 

It is hard to find the general form of distribution of cycle lengths for this situation.^ 
In the special case in which L is of the form 2", one can indeed see what will happen. 
For example, consider the structure shown in figure 11 if n = 4. Recall that the 
circuit structure makes the picture wrap around at j = 16, 32, ■ ■ ■ . After 32 steps, a 
structure arises between j = 32 and j = 48 which can cancel out the similar structure 
sitting between j = and j = 16. Then, after 48 steps, another structure beyond 
j = 48 unfolds itself and repeats the structure seen for small times. Hence the circuit 
has a closed cycle. The resulting picture is shown in Figure 15(B). Something very 
similar happens after 63 steps for L=7, but we have not fully unraveled the form of 
the resulting behavior for general L. 

The result shown in Table 3 partially answers a question which arose in our discussion 
of oscillations in the number of cycles for K = 2. In that case, we found that the 
preponderance of cycles had periods which were divisible by three. The problem is 
that neither linkages nor the trees put together from linkages have cycle lengths that 
are divisible by three. Large connected clusters might be expected to have a smoothly 
varying distribution of cycle lengths and hence only one third of these would be 
divisible by three. But small circuits can partially save the day. The table shows that 
the factor of three indeed arises quite often in these periods. We are pleased to see 
this, but not totally satisfied. There are only a few circuits (typically one) in a system 



Algebraic techniques introduced in reference [36] may be of use to solve the problem. 
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with A^ = 10 and K = 2. Somewhere there must be an additional source of factors 
of three. 
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5 Discussion 



In this paper we have studied the behavior of the Hamming distance and of the cycle 
length distribution in Boolean nets with time-reversible dynamics. We obtain strong 
evidence for the presence of a phase transition at a critical value K^ ^ 1.65 and 
present analytic bounds on K^. We observe large oscillations in the behavior of the 
Hamming distance as a function of time and of the cycle length distribution as a 
function of cycle length when K is relatively small. We propose that local structures, 
or small groups of spins with dynamics that do not depend on any other spins in the 
systems, play a crucial role in giving rise to these fluctuations. By analyzing local 
structures, significant analytic insight can be obtained regarding both the oscillatory 
behavior of the average number of cycle lengths averaged over realizations as well as 
the minima in the Hamming distance that occur at times that are multiples of 2" for 
integer n. 

Although we feel that we have a qualitative understanding of the phenomena exhib- 
ited by this model, some unexplored issues remain. It would be desirable to develop a 
more quantitative understanding of when period-three cycles are favored. Obtaining 
sharper bounds on the critical value Kc would also be desirable. 

We thank Raissa D'Souza for useful conversations. SNC and LPK gratefully acknowl- 
edge financial support from the National Science Foundation. SNC thanks the Aspen 
Center for Physics for hospitality during the preparation of this manuscript. 
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